###################################################
############ Supplemental Appendix Codes ##########
###################################################

rm(list=ls())

library(haven)
library(readr)
library(foreign)
library(ggplot2)
library(mvtnorm)
library(ggridges)

wd <- setwd("Your_Directory/JCR_Jooetal_replication")

############################
##### Figure A.2 - A.3 #####
############################

#Bootstrapped ATEs for India

indata <- read.csv("India_Data.csv")
set.seed(123456)
B <- 5000


#Function to Estimate Bootstrapped ATEs
bootTreat3 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=4)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    bootResults[i,1] <- mean(temp$EndorseForceQ11[which(temp$RWP==1)], na.rm=TRUE)
    bootResults[i,2] <- mean(temp$EndorseForceQ11[which(temp$Control==0)], na.rm=TRUE)
    bootResults[i,3] <- mean(temp$EndorseForceQ11[which(temp$Centrist==1)], na.rm=TRUE)
    bootResults[i,4] <- mean(temp$EndorseForceQ11[which(temp$Cooperative==1)], na.rm=TRUE)
    
    drop(list())
  }
  return(list(model=label, Control=bootResults[,2], RWP=bootResults[,1], Centrist=bootResults[,3], Cooperative=bootResults[,4]))
}

#Extract the Bootstrapped ATEs

##### Figure A2(a) #######
##### Full Sample #####
sepfull <- bootTreat3(indata, label="Full")

densT1 <- data.frame(y=sepfull$RWP, z=rep(c("Endorse Force"), each=B),  Group=rep(c("RW Populist"), each=B))
densT2 <- data.frame(y=sepfull$Centrist, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Centrist"), each=B))
densT3 <- data.frame(y=sepfull$Cooperative, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Cooperative"), each=B))

densC <- data.frame(y=sepfull$Control, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Control"), each=B))

dens1 <- rbind(densT1, densC)
dens2 <- rbind(densT2, densC)
dens3 <- rbind(densT3, densC)

dens_all<-rbind(densC,densT1,densT2,densT3)
dens_all$Group<-factor(dens_all$Group,levels = c("Control","RW Populist","Centrist","Cooperative"))

meansdens<-c(
  mean(densC$y),
  mean(densT1$y),
  mean(densT2$y),
  mean(densT3$y))

#Plot Density
pdf(paste(wd, "/Figures/figA2a.pdf", sep=""), width=10,height=6, paper='special')

ggplot(dens_all, aes(x = y, y = Group, fill = Group)) +
  geom_density_ridges() + xlab("")+ylab("")+
  theme_ridges() + 
  theme(legend.position = "none")+scale_fill_grey(start = .9, end = .7)+
  geom_point(aes(x=meansdens[1],y="Control"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[2],y="RW Populist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[3],y="Centrist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[4],y="Cooperative"),colour="black",size=3.5)+
  theme_minimal(base_size = 20)+ theme(legend.position = "none")

dev.off()


##### Figure A3(a) #######
####### India Partisan sample #######

indpart <- subset(indata, rpartisan==1)
seppart <- bootTreat3(indpart, label="Partisan")

densT1 <- data.frame(y=seppart$RWP, z=rep(c("Endorse Force"), each=B),  Group=rep(c("RW Populist"), each=B))
densT2 <- data.frame(y=seppart$Centrist, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Centrist"), each=B))
densT3 <- data.frame(y=seppart$Cooperative, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Cooperative"), each=B))

densC <- data.frame(y=seppart$Control, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Control"), each=B))

dens1 <- rbind(densT1, densC)
dens2 <- rbind(densT2, densC)
dens3 <- rbind(densT3, densC)

dens_all<-rbind(densC,densT1,densT2,densT3)
dens_all$Group<-factor(dens_all$Group,levels = c("Control","RW Populist","Centrist","Cooperative"))


meansdens<-c(
  mean(densC$y),
  mean(densT1$y),
  mean(densT2$y),
  mean(densT3$y))

pdf(paste(wd, "/Figures/figA3a.pdf", sep=""), width=10,height=6, paper='special')
ggplot(dens_all, aes(x = y, y = Group, fill = Group)) +
  geom_density_ridges() + xlab("")+ylab("")+
  theme_ridges() + 
  theme(legend.position = "none")+scale_fill_grey(start = .9, end = .7)+
  geom_point(aes(x=meansdens[1],y="Control"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[2],y="RW Populist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[3],y="Centrist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[4],y="Cooperative"),colour="black",size=3.5)+
  theme_minimal(base_size = 20)+ theme(legend.position = "none")
dev.off()

####### Japan Data #######

rm(list=ls())

wd <- setwd("Your_Directory/JCR_Jooetal_replication")

jpndata <-  read.csv("Japan_Data.csv")

set.seed(123456)
B <- 5000

bootTreat3 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=4)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    bootResults[i,1] <- mean(temp$EndorseForceQ11[which(temp$RWP==1)], na.rm=TRUE)
    bootResults[i,2] <- mean(temp$EndorseForceQ11[which(temp$Control==0)], na.rm=TRUE)
    bootResults[i,3] <- mean(temp$EndorseForceQ11[which(temp$Centrist==1)], na.rm=TRUE)
    bootResults[i,4] <- mean(temp$EndorseForceQ11[which(temp$Cooperative==1)], na.rm=TRUE)
    
    drop(list())
  }
  return(list(model=label, Control=bootResults[,2], RWP=bootResults[,1], Centrist=bootResults[,3], Cooperative=bootResults[,4]))
}

###### Figure A2(b) #######
##### Japan Full Sample #####
sepfull <- bootTreat3(jpndata, label="Full")

densT1 <- data.frame(y=sepfull$RWP, z=rep(c("Endorse Force"), each=B),  Group=rep(c("RW Populist"), each=B))
densT2 <- data.frame(y=sepfull$Centrist, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Centrist"), each=B))
densT3 <- data.frame(y=sepfull$Cooperative, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Cooperative"), each=B))

densC <- data.frame(y=sepfull$Control, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Control"), each=B))

dens1 <- rbind(densT1, densC)
dens2 <- rbind(densT2, densC)
dens3 <- rbind(densT3, densC)

dens_all<-rbind(densC,densT1,densT2,densT3)
dens_all$Group<-factor(dens_all$Group,levels = c("Control","RW Populist","Centrist","Cooperative"))


meansdens<-c(
  mean(densC$y),
  mean(densT1$y),
  mean(densT2$y),
  mean(densT3$y))


pdf(paste(wd, "/Figures/figA2b.pdf", sep=""), width=10,height=6, paper='special')
ggplot(dens_all, aes(x = y, y = Group, fill = Group)) +
  geom_density_ridges() +
  theme_ridges() + xlab("")+ylab("")+
  theme(legend.position = "none")+scale_fill_grey(start = .9, end = .7)+
  geom_point(aes(x=meansdens[1],y="Control"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[2],y="RW Populist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[3],y="Centrist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[4],y="Cooperative"),colour="black",size=3.5)+
  theme_minimal(base_size = 20)+ theme(legend.position = "none")
dev.off()

###### Figure A3(b) ########
####### Japan Partisan sample #######

jpnpart <- subset(jpndata, rpartisan==1)

seppart <- bootTreat3(jpnpart, label="Partisan")

densT1 <- data.frame(y=seppart$RWP, z=rep(c("Endorse Force"), each=B),  Group=rep(c("RW Populist"), each=B))
densT2 <- data.frame(y=seppart$Centrist, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Centrist"), each=B))
densT3 <- data.frame(y=seppart$Cooperative, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Cooperative"), each=B))

densC <- data.frame(y=seppart$Control, z=rep(c("Endorse Force"), each=B),  Group=rep(c("Control"), each=B))

dens1 <- rbind(densT1, densC)
dens2 <- rbind(densT2, densC)
dens3 <- rbind(densT3, densC)


dens_all<-rbind(densC,densT1,densT2,densT3)
dens_all$Group<-factor(dens_all$Group,levels = c("Control","RW Populist","Centrist","Cooperative"))

meansdens<-c(
  mean(densC$y),
  mean(densT1$y),
  mean(densT2$y),
  mean(densT3$y))


pdf(paste(wd, "/Figures/figA3b.pdf", sep=""), width=10,height=6, paper='special')
ggplot(dens_all, aes(x = y, y = Group, fill = Group)) +
  geom_density_ridges() + xlab("")+ylab("")+
  theme_ridges() + 
  theme(legend.position = "none")+scale_fill_grey(start = .9, end = .7)+
  geom_point(aes(x=meansdens[1],y="Control"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[2],y="RW Populist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[3],y="Centrist"),colour="black",size=3.5)+
  geom_point(aes(x=meansdens[4],y="Cooperative"),colour="black",size=3.5)+
  theme_minimal(base_size = 20)+ theme(legend.position = "none")
dev.off()


######### Parametric Model Figure Codes ##########

library(MASS)
library(mvtnorm)
library(ggplot2)
library(foreign)
library(dplyr)
library(ggplot2) 


#############################
########  Figure A4 #########
#############################

######## Figure A4(a) - India ########

rm(list=ls())

wd <- setwd("Your_Directory/JCR_Jooetal_replication")

idata <- read.csv("India_Data.csv")

men <- 1
age <- mean(idata$Age)
edu <- 5
marry <- 1
inc <- 3
auth <- 0
sup <- 1
trust <- 1
child <- 3
war <- 4
un <- 3
ideo <- mean(idata$IdeologyQ13)

B = 1000

coef <- c(1.315663,
          0.2331146,
          -0.1782619,
          -0.3789193,
          -0.0111952,
          0.0245004,
          0.2344052,
          -0.0918927,
          0.1901417,
          0.4327627,
          0.8503473,
          -0.1332827,
          -0.0649077,
          -0.1955854,
          0.2271978,
          0.005481
)

vcov <- read.csv("vcov_indfull.csv")

vcov <- as.matrix(vcov)
pr.bs1 <- matrix(NA,B,3)


set.seed(20221103)

for(i in 1:B){ 
  newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
  
  bs.sum0 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum1 <- newcoef[1]*1 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum2 <- newcoef[1]*0 + newcoef[2]*1 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum3 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*1 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  pr.bs1[i,1] <- plogis(bs.sum1) - plogis(bs.sum0)
  pr.bs1[i,2] <- plogis(bs.sum2) - plogis(bs.sum0)
  pr.bs1[i,3] <- plogis(bs.sum3) - plogis(bs.sum0)
  
}


nat <- data.frame(results=pr.bs1[,1], model="RW Populist")
m1 <- as.numeric(mean(nat$results))
u1 <- as.numeric(quantile(nat$results, p=0.995))
l1 <- as.numeric(quantile(nat$results, 0.005))

cent <- data.frame(results=pr.bs1[,2], model="Centrist")
m2 <- as.numeric(mean(cent$results))
u2 <- as.numeric(quantile(cent$results, 0.995))
l2 <- as.numeric(quantile(cent$results, 0.005))

coop <- data.frame(results=pr.bs1[,3], model="Cooperative")
m3 <- as.numeric(mean(coop$results))
u3 <- as.numeric(quantile(coop$results, 0.995))
l3 <- as.numeric(quantile(coop$results, 0.005))

CIdata <- data.frame(model=c("RW Populist", "Centrist", "Cooperative"), mean=c(m1, m2, m3), 
                     lower=c(l1, l2, l3), upper=c(u1,u2, u3))

par(mar=c(3,2,1,1))

pdf(paste(wd, "/Figures/figA4a.pdf", sep=""), width=5,height=4, paper='special')

ggplot() + 
  geom_errorbar(data=CIdata, mapping=aes(x=factor(model), ymin=lower, ymax=upper, group=factor(model)), width=0.2, size=.7, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=18),
                     axis.title.x = element_text(size=20),
                     axis.text.x = element_text(size=15),
                     axis.text.y = element_text(size=15))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.7) + 
  xlab("Treatment") + ylab("Change in Pr(Endorse Force)") + 
  scale_x_discrete(limits=c("RW Populist", "Centrist", "Cooperative"))

dev.off()


######## Figure A4(b) - Japan ########

rm(list=ls())

wd <- setwd("Your_Directory/JCR_Jooetal_replication")

jdata <- read.csv("Japan_Data.csv")

B <- 1000
men <- 1
age <- mean(jdata$Age)
edu <- 5
marry <- 1
inc <- 8
auth <- 0
sup <- 0
trust <- 1
child <- 0
war <- 3
un <- 2
ideo <- mean(jdata$IdeologyQ13)

coefjapan <- c(0.6929891,
               0.0247747,
               0.206073,
               0.0694105,
               -0.0012892,
               -0.0006472,
               0.2055089,
               0.0250849,
               0.1246571,
               0.2512075,
               0.229654,
               0.004229,
               0.2248228,
               0.0671294,
               0.1552802,
               -2.292278
)

vcovjapan <- read.csv("vcov_japanfull.csv")

vcovjapan <- as.matrix(vcovjapan)
pr.bs.japan <- matrix(NA,B,3)

B=1000
set.seed(20221103)

for(i in 1:B){ 
  newcoef <-  rmvnorm(1, mean=coefjapan, sigma=vcovjapan)
  
  bs.sum0 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum1 <- newcoef[1]*1 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum2 <- newcoef[1]*0 + newcoef[2]*1 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum3 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*1 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  pr.bs.japan[i,1] <- plogis(bs.sum1) - plogis(bs.sum0)
  pr.bs.japan[i,2] <- plogis(bs.sum2) - plogis(bs.sum0)
  pr.bs.japan[i,3] <- plogis(bs.sum3) - plogis(bs.sum0)
  
}


nat <- data.frame(results=pr.bs.japan[,1], model="RW Populist")
m1 <- as.numeric(mean(nat$results))
u1 <- as.numeric(quantile(nat$results, p=0.995))
l1 <- as.numeric(quantile(nat$results, 0.005))

cent <- data.frame(results=pr.bs.japan[,2], model="Centrist")
m2 <- as.numeric(mean(cent$results))
u2 <- as.numeric(quantile(cent$results, 0.995))
l2 <- as.numeric(quantile(cent$results, 0.005))

coop <- data.frame(results=pr.bs.japan[,3], model="Cooperative")
m3 <- as.numeric(mean(coop$results))
u3 <- as.numeric(quantile(coop$results, 0.995))
l3 <- as.numeric(quantile(coop$results, 0.005))

CIdata <- data.frame(model=c("RW Populist", "Centrist", "Cooperative"), mean=c(m1, m2, m3), 
                     lower=c(l1, l2, l3), upper=c(u1,u2, u3))

pdf(paste(wd, "/Figures/figA4b.pdf", sep=""), width=5,height=4, paper='special')

par(mar=c(3,2,1,1))

ggplot() + 
  geom_errorbar(data=CIdata, mapping=aes(x=model, ymin=lower, ymax=upper), width=0.2, size=.7, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=18),
                     axis.title.x = element_text(size=20),
                     axis.text.x = element_text(size=15),
                     axis.text.y = element_text(size=15))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), linewidth=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.7) + 
  xlab("Treatment") + ylab("Change in Pr(Endorse Force)")+
  scale_x_discrete(limits=c("RW Populist", "Centrist", "Cooperative"))

dev.off()



#############################
######## FIGURE A5 ##########
#############################

#### Figure A5(a) - India R-wing Partisan ####

rm(list=ls())

wd <- setwd("Your_Directory/JCR_Jooetal_replication")

idata <- read.csv("India_Data.csv")

men <- 1
age <- mean(idata$Age)
edu <- 5
marry <- 1
inc <- 3
auth <- 0
sup <- 1
trust <- 1
child <- 3
war <- 4
un <- 3
ideo <- mean(idata$IdeologyQ13)

B = 1000

coef1 <- c(1.141048,
           0.6475108,
           0.0748971,
           -0.1695816,
           -0.0374145,
           -0.3195208,
           -0.128713,
           -0.1458982,
           0.0405899,
           0.9226759,
           1.758717,
           -0.0980382,
           -0.3139624,
           0.2463492,
           0.4505681,
           -0.1033389
)

vcov1 <- read.csv("vcov_indiarpart.csv")

vcov1 <- as.matrix(vcov1)

pr.bs2 <- matrix(NA,B,3)


set.seed(20221103)

for(i in 1:B){ 
  newcoef <-  rmvnorm(1, mean=coef1, sigma=vcov1)
  
  bs.sum0 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum1 <- newcoef[1]*1 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum2 <- newcoef[1]*0 + newcoef[2]*1 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum3 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*1 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  pr.bs2[i,1] <- plogis(bs.sum1) - plogis(bs.sum0)
  pr.bs2[i,2] <- plogis(bs.sum2) - plogis(bs.sum0)
  pr.bs2[i,3] <- plogis(bs.sum3) - plogis(bs.sum0)
  
}


nat1 <- data.frame(results=pr.bs2[,1], model="RW Populist")
m11 <- as.numeric(mean(nat1$results))
u11 <- as.numeric(quantile(nat1$results, p=0.995))
l11 <- as.numeric(quantile(nat1$results, 0.005))

cent1 <- data.frame(results=pr.bs2[,2], model="Centrist")
m21 <- as.numeric(mean(cent1$results))
u21 <- as.numeric(quantile(cent1$results, 0.995))
l21 <- as.numeric(quantile(cent1$results, 0.005))

coop1 <- data.frame(results=pr.bs2[,3], model="Cooperative")
m31 <- as.numeric(mean(coop1$results))
u31 <- as.numeric(quantile(coop1$results, 0.995))
l31 <- as.numeric(quantile(coop1$results, 0.005))

CIdata1 <- data.frame(model=c("RW Populist", "Centrist", "Cooperative"), mean=c(m11, m21, m31), 
                      lower=c(l11, l21, l31), upper=c(u11,u21, u31))

par(mar=c(3,2,1,1))

pdf(paste(wd, "/Figures/figA5a.pdf", sep=""), width=5,height=4,paper='special')

ggplot() + 
  geom_errorbar(data=CIdata1, mapping=aes(x=model, ymin=lower, ymax=upper, group=factor(model)), width=0.2, size=.7, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=18),
                     axis.title.x = element_text(size=20),
                     axis.text.x = element_text(size=15),
                     axis.text.y = element_text(size=15))  +
  geom_point(data=CIdata1, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.7) + 
  xlab("Treatment") + ylab("Change in Pr(Endorse Force)") +
  scale_x_discrete(limits=c("RW Populist", "Centrist", "Cooperative"))


dev.off()


#### Figure A5(b) - Japan R-wing Partisan ####

rm(list=ls())

wd <- setwd("Your_Directory/JCR_Jooetal_replication")

jdata <- read.csv("Japan_Data.csv")

B <- 1000
men <- 1
age <- mean(jdata$Age)
edu <- 5
marry <- 1
inc <- 8
auth <- 0
sup <- 0
trust <- 1
child <- 0
war <- 3
un <- 2
ideo <- mean(jdata$IdeologyQ13)

coefjapan1 <- c(2.409163,
                0.9096439,
                1.151435,
                0.4855414,
                0.0655983,
                0.353038,
                -1.661666,
                0.1298928,
                -0.0758032,
                0.1837595,
                1.719275,
                0.1222071,
                1.095577,
                0.3804429,
                0.0979223,
                -11.66576)

vcovjapan1 <- read.csv("vcov_japanrwpart.csv")

vcovjapan1 <- as.matrix(vcovjapan1)

pr.bsjapan2 <- matrix(NA,B,3)


set.seed(20221103)

for(i in 1:B){ 
  newcoef <-  rmvnorm(1, mean=coefjapan1, sigma=vcovjapan1)
  
  bs.sum0 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum1 <- newcoef[1]*1 + newcoef[2]*0 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum2 <- newcoef[1]*0 + newcoef[2]*1 + newcoef[3]*0 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  bs.sum3 <- newcoef[1]*0 + newcoef[2]*0 + newcoef[3]*1 + newcoef[4]*men + newcoef[5]*age +newcoef[6]*edu + 
    newcoef[7]*marry + newcoef[8]*inc + newcoef[9]*auth + newcoef[10]*sup + newcoef[11]*trust + newcoef[12]*child + 
    newcoef[13]*war + newcoef[14]*un + newcoef[15]*ideo + newcoef[16]
  
  pr.bsjapan2[i,1] <- plogis(bs.sum1) - plogis(bs.sum0)
  pr.bsjapan2[i,2] <- plogis(bs.sum2) - plogis(bs.sum0)
  pr.bsjapan2[i,3] <- plogis(bs.sum3) - plogis(bs.sum0)
  
}


nat1 <- data.frame(results=pr.bsjapan2[,1], model="RW Populist")
m11 <- as.numeric(mean(nat1$results))
u11 <- as.numeric(quantile(nat1$results, p=0.995))
l11 <- as.numeric(quantile(nat1$results, 0.005))

cent1 <- data.frame(results=pr.bsjapan2[,2], model="Centrist")
m21 <- as.numeric(mean(cent1$results))
u21 <- as.numeric(quantile(cent1$results, 0.995))
l21 <- as.numeric(quantile(cent1$results, 0.005))

coop1 <- data.frame(results=pr.bsjapan2[,3], model="Cooperative")
m31 <- as.numeric(mean(coop1$results))
u31 <- as.numeric(quantile(coop1$results, 0.995))
l31 <- as.numeric(quantile(coop1$results, 0.005))

CIdata1 <- data.frame(model=c("RW Populist", "Centrist", "Cooperative"), mean=c(m11, m21, m31), 
                      lower=c(l11, l21, l31), upper=c(u11,u21, u31))

pdf(paste(wd, "/Figures/figA5b.pdf", sep=""), width=5,height=4,paper='special')

par(mar=c(3,2,1,1))

ggplot() + 
  geom_errorbar(data=CIdata1, mapping=aes(x=model, ymin=lower, ymax=upper), width=0.2, size=.7, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=18),
                     axis.title.x = element_text(size=20),
                     axis.text.x = element_text(size=15),
                     axis.text.y = element_text(size=15))  +
  geom_point(data=CIdata1, mapping=aes(x=model, y=mean), size=2,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.7) + 
  xlab("Treatment") + ylab("Change in Pr(Endorse Force)") +
  scale_x_discrete(limits=c("RW Populist", "Centrist", "Cooperative"))

dev.off()


########################
###### Figure A.7 ######
########################

##### Figure A.7(a) #####

mean=c(.1193803, .0423279)  # Estimates based on the margins code in STATA
se=c(.0585588 ,  .033742)  # Estimates based on the margins code in STATA

CIdata <- data.frame(model=c("No Rival", "Rival"), mean=mean, 
                     lower=c(mean[1]-1.96*se[1],mean[2]-1.96*se[2] ), 
                     upper=c(mean[1]+1.96*se[1],mean[2]+1.96*se[2] ))

ggplot() + 
  geom_errorbar(data=CIdata, mapping=aes(x=model, ymin=lower, ymax=upper), width=0.2, size=1, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=14, face="bold"),
                     axis.text.x = element_text(size=14, face='bold'),
                     axis.text.y = element_text(size=10))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2.5,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.7) + 
  xlab("") + ylab("Marginal Effect of RWP") 


CIdata <- data.frame(model=c("No Rival", "Rival"), mean=mean, 
                     lower=c(mean[1]-2.576*se[1],mean[2]-2.576*se[2] ), 
                     upper=c(mean[1]+2.576*se[1],mean[2]+2.576*se[2] ))

pdf(paste(wd, "/Figures/figA7a.pdf", sep=""), width=6,height=6,paper='special')

ggplot() + 
  geom_errorbar(data=CIdata, mapping=aes(x=model, ymin=lower, ymax=upper), width=0.2, size=1, color="black") + 
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
                     axis.title.y = element_text(size=14, face="bold"),
                     axis.text.x = element_text(size=14, face='bold'),
                     axis.text.y = element_text(size=10))  +
  geom_point(data=CIdata, mapping=aes(x=model, y=mean), size=2.5,  fill="white") +
  geom_hline(yintercept=0, linetype="dashed",  size=0.7) + 
  xlab("") + ylab("Marginal Effect of RWP") 

dev.off()


##### Figure A.7(b) #####
wd <- setwd("Your_Directory/JCR_Jooetal_replication")

coef <- c(2.877169,
          2.103525,
          -1.172508,
          -0.1707873,
          -0.1595341,
          -0.7620533,
          4.61949,
          5.968458,
          0.9610305,
          0.8595328,
          0.4809819,
          0.522015,
          0.0624043,
          0.2383076,
          -0.379984,
          -0.0509925,
          0.1392838,
          0.0719939,
          1.350758,
          -10.45478)

vcov <- as.matrix(read.csv("vce_rnrdyadicdum.csv"))

dat<- read.dta("dyadicdata.dta")


lag <- 0
c1 <- median(dat$cinc1, na.rm = T)
c2 <- median(dat$cinc2, na.rm = T)
p2 <- median(dat$partipdem2, na.rm=T)
r2 <- 0 
majp1 <- 0
majp2 <- 0 
ally <- 0
male1 <- 1
edu1 <- 3
mil1 <- 0
milcar <- 0
combat <- 0
reb1 <- 0
ireg1 <- 0
rival <- 0


B <- 1000

pardem <- seq(0,1, by =0.01) 

out1 <- matrix(NA, ncol=3, nrow=length(pardem))

# pinit1_lag	rightparti1	partipdem1	partipdem2	right1	right2	cinc1	cinc2	majpower1	majpower2	
# atopally	male1	leveledu1	milcombat1	militaryca~1	combat1	rebel1	irregular1	rival_dum	_cons


for(i in 1:length(pardem)){ 
  pdem <- pardem[i] 	
  pr.bs <- matrix(NA,B,1)
  for(b in 1:B){ 
    newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
    bs.sum <- newcoef[1]*lag + newcoef[2]*0*pdem + newcoef[3]*pdem + newcoef[4]*p2 + newcoef[5]*0 +newcoef[6]*r2 + newcoef[7]*c1
    + newcoef[8]*c2 + newcoef[9]*majp1 + newcoef[10]*majp2 + newcoef[11]*ally + newcoef[12]*male1 + newcoef[13]*edu1
    + newcoef[14]*mil1 + newcoef[15]*milcar + newcoef[16]*combat + newcoef[17]*reb1 + newcoef[18]*ireg1 + newcoef[19]*rival + newcoef[20]
    
    
    bs.sum2 <-  newcoef[1]*lag + newcoef[2]*1*pdem + newcoef[3]*pdem + newcoef[4]*p2 + newcoef[5]*1 +newcoef[6]*r2 + newcoef[7]*c1
    + newcoef[8]*c2 + newcoef[9]*majp1 + newcoef[10]*majp2 + newcoef[11]*ally + newcoef[12]*male1 + newcoef[13]*edu1
    + newcoef[14]*mil1 + newcoef[15]*milcar + newcoef[16]*combat + newcoef[17]*reb1 + newcoef[18]*ireg1 + newcoef[19]*rival + newcoef[20]
    
    pr.bs[b,] <- plogis(bs.sum2) - plogis(bs.sum) 
  }
  out1[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE)) 
  rm(bs.sum)
}

pdf(paste(wd, "/Figures/figA7b.pdf", sep=""), width=6,height=6,paper='special')

par(mar=c(5,5,1.5,1.5))
plot(pardem, out1[,1], type="n", xlab="Participatory Democracy", ylab="Change in Pr(MIDS=1)",
     cex.lab=1.5, cex.axis=1.5, ylim = c(-0.3,1))
lines(lowess(pardem, out1[,1]), lty=1, lwd=2)
lines(lowess(pardem, out1[,2]), lty=2, lwd=2)
lines(lowess(pardem, out1[,3]), lty=2, lwd=2)
abline(h=0, lty=2, col='red', lwd=1)

dev.off()


# End of File 